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Nonlinear evolution of three-dimensional strained islands or quantum dots in heteroepitaxial thin 
films is studied via a continuum elasticity model and the development of a nonlinear dynamic 
equation governing the film morphological profile. All three regimes of island array evolution are 
identified and examined, including a film instability regime at early stage, a nonlinear coarsening 
regime at intermediate times, and the crossover to a saturated asymptotic state, with detailed 
behavior depending on film-substrate misfit strains but not qualitatively on finite system sizes. The 
phenomenon of island stabilization and saturation, which corresponds to the formation of steady 
CU but non-ordered arrays of strained quantum dots, occurs at later time for smaller misfit strain. It is 

found to be controlled by the strength of film-substrate wetting interaction which would constrain 
the valley-to-peak mass transport and hence the growth of island height, and also determined by the 
, effect of elastic interaction between surface islands and the high-order strain energy of individual 

^s^j ■ islands at late evolution stage. The results are compared to previous experimental and theoretical 

studies on quantum dots coarsening and saturation. 
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I. INTRODUCTION 



The formation of surface nanostructures such as islands or quantum dots during heteroepitaxy of strained films 
has attracted continuing great interest, due to its importance in both fundamental understanding of material growth 
and its wide range of applications particularly for optoelectronic nano devices^—. One of the underlying mechanisms 
of island formation has been attributed to the occurrence of morphological instability of the strained film (i.e., the 
Asaro-Tiller-Grinfeld instability^), for which the competition between the stabilization effect of film surface energy 
and the destabilization effect of system elastic energy due to film-substrate misfit strain plays a key role. This 
involves the procedure of stress release in the film via surface mass transport driven by local energy gradient and 
consequently the formation of coherent nonplanar surface structures like nanoscale islands^. Compared to other 
formation mechanisms such as thermally activated nucleation, a main feature here is the continuous, nucleationless 
nature of the film instability and the subsequent island growth that starts from rough surface 7 - - — . It results in more 
regular and correlated arrays of self assembled quantum dots that are appealing for a variety of applications. 
. The onset of this type of quantum dots or strained islands formation has been well understood, as studied via detailed 
1 instability analysis for both single-component and alloy strained films 10 . - .! 3 , as well as multilayers /superlattice o 14 i 1 5 . 
t-H . Such film instability, showing as surface ripples and cell-like undulations which are the precursor of coherent island 
formation, has been observed in experiments of semiconductor heteroepitaxial films such as SiGe 7 -^. On the other 
hand, the subsequent nonlinear evolution of these strained dots that are well beyond the initial linear stability stage, 
H 1 is more complicated and much less understood. A typical phenomenon is the coarsening of quantum dot islands, 
r—i | showing as the increase of average island size during film evolution and the shrinking of small dots. Such scenario 
has been observed in experiments of Ge/Si(001)i^, SiGe/Si(001)22, and InAs/GaAs(001) 21 , although with different 
mechanisms and behavior of coarsening reported. For the example of Ge-Si systems that have been extensively studied, 
nonlinear island coarsening rate that deviates from classical results of Ostwald ripening has been found; however, it 
was associated with different mechanisms of island morphology evolution in different experiments, such as the effect 
of island shape transition accompanying the coarsening process^— or a contrasting kinetic picture that incorporates 
clastic interactions but not island shape transition^. Also, the slowing^ or suppression^ of coarsening process at 
late stage and the resulting saturated, stabilized quantum dot array a 17 i 19 ' 22 have been observed in some experiments, 
but not others 2 ^. 

The corresponding theoretical/computational study on nonlinear island evolution is also far from conclusive. Most 
efforts are based on continuum approaches including continuum elasticity theory 2 ^— and phase field method o 37 ' 38 , 
in addition to other modeling techniques that incorporate crystalline details, such as kinetic Monte Carlo method 
with elastic interaction 3 ^— and the recently developed phase-field-crystal (PFC) model and the associated amplitude 
equation formalism^ - —. For continuum elasticity modeling, which is the current main avenue for studying strained 
island coarsening due to the large length and time scales involved, much recent focus has been put on the derivation 
and simulation of nonlinear evolution equations through approximating the system elasticity and the dynamics of film 
morphology via perturbation methods. Two limits of system configuration were addressed in early studies of such 
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approach, including the limit of perfectly rigid substrat o 24 ' 30 and the case of infinitely thick strained film 2 ^. More 
recently, similar approximation has been applied to heteroepitaxial systems consisting of a strained thin film grown on 
an elastic substrate as configured in most experiments. The corresponding reduced nonlinear evolution equations have 
been derived and simulated, based on the long-wave or small-slope approximation of film surface profil o 31 ' 34 ' 35 or the 
assumption of small surface gradient 33 . Some physical mechanisms in thin film growth have also been incorporated, 
such as the wetting effect between the strained film and the underlying substrate. These approaches are more efficient 
for large scale simulations, as compared to directly solving the full system elasticity problem and the corresponding 
full dynamic equation of morphological profile, which instead is of high computational cost and hence usually involves 
limited system size, island number, and evolution time particularly for three-dimensional (3D) systems. 

Despite the success of these theoretical approaches in describing properties of quantum dot formation and film 
morphology, some behavior of film nonlinear evolution, particularly the process of island coarsening vs. saturation, is 
still not well understood, with inconsistent results given in different studie s 26 ! 29 ' 32- — . For the case of film annealing 
as examined in most simulations, coarsening of strained island arrays has been reproduced, although with different 
coarsening rate found in different approache d 29 ! 34 ' 35 . One of the main difference in these work is the result for 
asymptotic and steady state of the film morphology. Stable arrays of quantum dots that persist after the coarsening 
stage have been obtained in both studies of the reduced nonlinear evolution equatio n 33 ' 36 and the direct solution of 
the full elasticity proble m 26 ' 32 , consistent with the observation in some Si-Ge experiment o 17 ' 19 ' 22 . However, as in some 
other experiments 2 - such scenario of the suppression or cessation of island coarsening was not found in other modeling 
processe s 29 ' 34 , and recent nonlinear analysis of a system evolution equation suggested that a regular quantum dot 
array would be unstable as a result of subcritical bifurcation^ 4 -. Note that in many previous studies the saturation 
of coarsening islands has been attributed to the effect of surface energy anisotrop y 26 ' 32 ' 36 , although the saturation 
phenomenon has also been observed in recent simulation without such anisotropy effect 3 - 3 ". These discrepancies in 
modeling results could probably be related to different types of approximation and various ways of small variable 
expansion and truncation involved in the approaches, and/or the difficulty in simulating large enough system size and 
long enough evolution time required for experimental comparison. 

In this paper we focus on the nonlinear evolution of strained quantum dot islands grown epitaxially on an elastic 
substrate, based on a continuum elasticity model and the development of a systematic approach for approximately 
solving the film-substrate elastic state via a perturbation analysis in Fourier space. Results up to second- and third- 
order perturbation of surface morphology are presented, and our approach can be readily extended to incorporate 
higher-order solutions. We can then derive a new nonlinear evolution equation governing the dynamics of strained film 
morphology, which allows us to systematically examine the detailed behavior of island evolution at large enough spatial 
and temporal scales. This nonlinear equation, with the incorporation of the wetting effect and also a second-order 
truncation in the elastic solution, is applied to the study of the coarsening and saturating process in post-deposited, 
annealing films. Our focus is on systems of small misfit strains, which correspond to large enough length scale of the 
resulting surface structure as compared to the scale of crystalline lattice. (This is based on recent studies^ 3 - showing 
that continuum approaches, such as the continuum elasticity theory developed here, can well describe the films in 
weak strain limit, but not for large misfit stress which would lead to qualitatively different behavior of surface islands 
due to the effects of discrete lattice structure.) The whole series of island evolution can be reproduced in our numerical 
simulations, including three characteristic regimes: the development of morphological instability and island formation 
at early times, nonlinear coarsening of islands at intermediate stage, and the slowing of such coarsening process 
which leads to a saturated state of steady quantum dot arrays. We also investigate the mechanisms underlying the 
phenomenon of island stabilization and saturation, based on the study of wetting effect on the constraint of surface 
mass transport between island valleys and peaks, and also on a detailed examination of strain relaxation process via 
studying the temporal evolution and spatial distribution of elastic energy density at the film surface. We identify a 
new factor responsible for suppressing the island growth and coarsening, which is attributed to the effect of high-order 
elastic energy of individual islands and the elastic interaction between them. 

II. MODEL 

Assume that a strained film of spatially varying height h(x,y,t) is deposited epitaxially on a semi-infinite elastic 
substrate that occupies the region z < 0. The misfit strain in the film with respect to the substrate is given by 
e = (a; — a s )/a s , where a/ and a s are the lattice spacings of the film and the substrate respectively. For such 
coherent, dislocation-free system, the evolution of the film surface morphological profile h(x,y,t) is governed by 



dh 2 5T 



(1) 
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where Th is the kinetic coefficient determined by surface diffusion, V 2 is the surface Laplacian, v is the deposition 
rate, and g = 1 + |V/i| 2 . Here the effect of film-substrate interdiffusion is neglected. The total free energy functional 

J- consists of two parts, including the elastic energy T e i — f_ d 3 r£ where £ represents the strain energy density, 
and the surface free energy T s = J d 2 rj s (h)y/g where 7 S is the thickness-dependent, isotropic surface tension with 
the effect of wetting interaction between the film and substrate incorporated. The dynamic equation ([1]) can then be 
rewritten as£ ' 31 i 34 



= T h JgV 2 s [ 7 /t + W{h) + £f] + v, 



(2) 



where k is the mean surface curvature, 7 is the surface energy density, £ * gives the film elastic energy density at the 
surface z = h, and the wetting potential W can be approximated via a phenomenological glued-layer wetting model 34 
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Here w gives the strength of the film-substrate wetting interaction, h m i is the characteristic wetting-layer thickness 
that is usually of few monolayers, and the exponent a w (> 0) gives the singularity of the potential W in the limit of 
h — > that emulates the persistence of the wetting layer during film evolution 3 ^. 

The formulation of elasticity for this film-substrate system has been well established^. In isotropic, linear elasticity 
theory (i.e., with harmonic approximation), the elastic energy density is given by £ = haijUij, where i,j = x,y,z, 



and Uij is the linear elastic strain tensor defined by it^ = 
Hooke's law for isotropic elastic system, the stress tensor cr, 



(djUi + diUj)/2 (with m the displacement field). 
j in the strained film is expressed by 
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where fi is the shear modulus and v is the Poisson ratio. The stress tensor in the substrate is also given by Eq. Q 
with e = 0. Here for simplicity we have assumed equal elastic constants in the film and substrate, which corresponds 
to the situation in most experimental systems where the difference of elastic constants between film and substrate 
materials is not significant. 

Since the elastic relaxation occurs at a time scale of orders of magnitude faster than that of the atomic diffusion 
process and the associated system morphological evolution, it is usually assumed that the mechanical equilibrium 
condition djdij — is always satisfied in both film and substrate. Using Eq. (j4} we can obtain Navier's equations in 
the whole film-substrate system 

(f - 2v)d 2 jUl + didjUj = 0. (5) 
The corresponding boundary condition on the top film surface is given by 

at z = h(x, y, t), (6) 



due to the neglecting of external pressure on the free surface, while the boundary conditions at the film-substrate 
interface is determined by the continuity of stress and displacement fields: 



at z 



0. 



(7) 



Here rij is the unit vector normal to the film surface or the film-substrate interface, and the subscripts "/" and "s" 
refer to the film and substrate phases, respectively. Also, inside the substrate region which is far away from the film 
we have 



-> 



for 



(8) 



III. PERTURBATION ANALYSIS AND NONLINEAR EVOLUTION EQUATION 



To solve this elasticity problem, we adopt a perturbation analysis in Fourier space based on the expansion of the 
small vertical variation of film surface profile. More specifically, given the Fourier transform of the film morphological 
profile 



h = h + y]h(q, t)( 



(9) 
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where h = ho + vt is the average film thickness at any time t (with ho the initial film thickness), the Fourier components 
of the displacement field Uj(q), stress tensor Oij(q) (i, j = x, y, z), and film elastic energy density £■> are expanded in 
the order of surface perturbation h(q), i.e., 



Ui = Ui + 



,i(q x x+q y y) 



.(1) . .(2) 

Ui = u. + u) 



-(3) 



£ s = £ f +J2 £ f (q)e l ^ x +^ 



. .(2) . »(3) 



(10) 



y) 



f _ 



For the Oth-order base state with planar, uniformly strained film, the elasticity solutions are given by£: u 
(i = x,y) and u{ — u{ z z for the displacement fields, the strain tensor u{j — except for u{ z = e(l +u)/(l — v), the 
stress tensor a(- = except for &l x = a^y = —2/j,u® z , and the Oth-order elastic energy density — Ee 2 v) (where 
E is the Young's modulus). For the substrate, the corresponding base state is stress- free, with u\ = ufj = afj = 0. 

The elastic properties at higher orders can be obtained by substituting the expansions © and (|10[) into Eqs. 
In Fourier space the Navier's equations ([3]) can be rewritten as 



(1 - 2v){dl - q 2 )uf + i Qj [iq x U& + iq v U^ + d z &® 
(1 - 2v){d 2 z - q 2 )u { P + d z \iq x u^ + iq yU f + d z u { P 



= 0, 

= o, 



for j = x, y, 



(11) 
(12) 



for £th order expansion (£ = 1, 2, 3, ...). The corresponding general solutions have the same format as that obtained 
in Ref. 6 for 1st order equations, which read 
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for the film, and 
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(14) 



for the substrate after using the boundary conditions (J7J and (j5J) at the film-substrate interface and inside the 
substrate. The coefficients af\ (3^ ] , C®, D®, and in Eqs. CE3]) and (JHJ) are determined via the expansion of 
boundary conditions (©-([7]) in orders of perturbation h. Note that the 1st order solution has been known with the 
use of linearized boundary condition o 6 ' 10 , with the perturbed elastic energy density being given by 



£ (Df = _lE{l + v) ^ 2 



e 2 qh(q). 



(15) 



For the 2nd order expansion of the boundary conditions, at the top surface of the film, z = h, Eq. ([6]) can be 
written as 

- - O^ff' (q')fc(q -q')-E - O^V^q - q') + ^^(q) = °> 

q' q' 



(16) 



while the continuity of stress at the film-substrate interface z — (i.e., Eq. (O) yields 

-(2)// \ -(2)s/ \ 

°)z (q) = o-}» (q). 

with j = x,y,z. Substituting Eqs. (IIUI) and (ITU to these boundary conditions (IT^l) and (IT71) . the second order 
coefficients of the solution can be obtained as follows: 
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where 



and 



QxQ' 



+2Ee^ [q x (q y - q' y ) + q y (q x - q' x )} 
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(19) 



(20) 



Based on the above solution, we can determine the second order elastic energy density, i.e. 
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where the expressions of lst-order stress tensor fly at the top surface are given in the appendix [see Eqs. (IA.11|) 
(|A.16I) ]. The 2nd-order stress tensor can be also calculated, with results shown in Eqs. (|A.5|) — (|A.10I) of the appendix. 

We then derive the nonlinear evolution equation for film surface morphology from Eq. @, using the results of 
perturbed elastic energy density given in Eqs. (|15[) and (|2ip . All the terms in the dynamic equation @ are expanded 
up to second order of surface perturbation h, except for W(h) for which the full nonlinear wetting potential form 
Eq. ([3]) is used. For the surface energy term jk, noting that the surface curvature k = — V • [V/i/yl + | V/i| 2 ], 
we have 7K ~ —j'S/ 2 h + 0(h 3 ) and hence only need to keep the linear order term in the 2nd-order approximation 
considered here. To further simplify the calculation, we choose a length scale I = 7/^0 and a time scale t — l A /'yTh, 
where the characteristic strain energy density £q — 2Ee^{\ + v)j(X — v) with eo a reference misfit value. The resulting 
nondimensional dynamic equation for the perturbed surface profile h(q,t) is given by 



dh 



(-<Z 4 + e*Y)ft 



q 2 W q 
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q')A(q,q'), 



(22) 



where W q is the Fourier transform of the rescaled wetting potential W(h)/£o, e* = c/cq, and 



A(q,q') =q 2 
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q'l 
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In Eq. (|22[) . the first term of the right-hand-side is the combination of the surface energy contribution and the 1st order 
clastic energy density £ 'W, consistent with the previous linear-order results ^ 31 ' 34 . The last term is from the 2nd-order 
perturbation result of the elastic energy density, i.e., Eq. (f2"Tj) for £ ^>f . It would be straightforward, although with 
more complicated processes, to extend the above approach to incorporate higher-order contributions, based on the 
perturbed analysis of system elasticity given in Eqs. (|9))- (|14[) . That is, to obtain the nth-order elastic results, we can 
first express the nth-order expansion of the boundary conditions in terms of £th-order (£ = l,2,...,n — l,n) stress 
tensors [similar to the expression in Eq. (|16l) ]. and use it to calculate the nth-order solution of the displacement field 
given in Eq. (|13[) ; The corresponding film clastic properties can then be derived, particularly the nth-order elastic 
energy density £( n >f at the film surface which can be expressed as a function of £th-order (£ = 1, 2, n — 1) elastic 
quantities that are already known. We have applied this recursive method to third-order calculations, with results 
of elastic energy density £( 3 )^ shown in the appendix. The corresponding higher-order evolution equation can be 
obtained via adding term —q 2 £^ 3 ^ to Eq. (|22|) . 



IV. LINEAR STABILITY ANALYSIS 



In the following studies of film evolution and strained island dynamics, we focus on Eq. (|2"2")l with 2nd-order elastic 
properties. Here we first perform a linear stability analysis of Eq. (|22[) to determine the conditions of morphological 
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instability of the system; such conditions are needed for the nonlinear calculations given in Sec. |V] and Sec. I VII 
Following the standard procedure, we assume an exponential growth h = /loexp(cr^i) at early time, and apply it to 
the linearized evolution equation of h. The characteristic equation for the perturbation growth rate ah is then given 
by 



Oh 



-g 4 + gV 2 



g 2 a, 



(24) 

h m i/l. From the above dispersion 
relation we can identify the condition for the occurrence of film morphological instability, which is given by e* > 4a 
or equivalently, 



where a — (w* /h*,)(x + a w )x Q,u 1 e x , x — h/h m i, w* = w/£q, and h 
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The corresponding characteristic wave number of film instability (for the fastest instability growth mode) can be 
written as 



*2 



:-S2 



(26) 



Eq. (|25[) is used to identify the parameters in our numerical simulations shown below, for which the initial film 
instability and hence the appearance of nonplanar surface morphology or islands are required. Note that for given 
film conditions such as misfit strain e and wetting parameters, Eq. (1251) indicates that due to the film-substrate 
wetting effect (with a w > 0), the morphological instability and surface nanostructures will develop only for thick 



enough films, with the critical thickness h c determined by Eq. (|25p. i.e., e x 



/(X C 



Aw* /e* h* ml (where 



x c = h c /h m i); the value of h c increases with smaller film-substrate misfit strain. Also, the characteristic size (or 
wavelength A = 27r/g max ) of surface structures at the initial stage will decrease with the increasing average film 
thickness h, as can be obtained from Eq. (|26|) . 



V. EFFECTS OF WETTING POTENTIAL 



To validate our model system and the nonlinear dynamic equation (|22[) derived above, we first examine the effect of 
wetting potential on film evolution and compare it to the well-known results of cusp formation obtained from previous 
full elasticity calculation ;] 23 ' 25 . Since Eq. (|22l) is already presented in Fourier space, in our numerical simulations we 
directly use the spectral method with periodic boundary conditions along the lateral x and y directions, and also 
an exponential propagation algorithm for time integration 46 . This allows us to use large enough time steps (up to 
At = 1 in most results shown below, except for the study of groove/cusp formation for which At = 0.01 is used). 
For rescaling parameters that are associated with the nondimensional equation (|22j) . the reference misfit eo = 3% is 
chosen, and thus the length scale can be estimated as I ~ 5.5 nm if using the material parameters of Ge/Si system. 
For simplicity, in this work we focus on the case of nongrowing films with deposition rate v = and simulate the 
annealing process of film evolution. 

As found first by Yang and Srolovitz 2 ^, and later in various numerical studies of 2D 25 ' 28 and 3D 30 ' 33 systems via 
solving either the full elasticity problem or the reduced nonlinear evolution equations, deep grooves or cusps will form 
in stressed solid systems without the incorporation of wetting effect. This is well reproduced in our numerical results 
of Eq. (|22|) . as shown in Fig. [T]for 3% misfit films with v = 1/3 and initial film thickness ho = 0. Two types of 
initial conditions are used: (1) a small random disturbance of a planar film of thickness ho, with results of a 128 x 128 
system presented in Fig. d](a) and (b), and (2) a weakly perturbed film with doubly-periodic sinusoidal surface profile 
h = ho + Ao[cos(q x ax) + cos(q y oy)], as given in panels (c) and (d) of Fig. Q]for a A^o X A y o system (where A^o = 2ir/q x o 
and X y o — 2ir/q y o). For condition (2) the initial perturbed amplitude is set as Aq = 0.01, and a perturbed wavevector 
1x0 = QyO = 3/4y2 is chosen, corresponding to the wave number of the most linearly unstable mode determined by 
Eq. (|26[) . In both cases the formation of singular cusps or deep grooves and their rapid growth have been found 
during the film evolution, as evidenced by the 3D morphological profiles given in Fig. Q] (a) and (c), and also from 
the results of time-evolving 2D cross-section profiles shown in (b) and (d) which are consistent with the previous 2D 
results of Spencer and Meiror*^ 5 . and Xiang and Fj2£. 

To incorporate the wetting effect, in our calculations we use a pseudospectral method; that is, we first evaluate 
the wetting potential W(h) from Eq. ((3]) in real space and then obtain its Fourier component W q as used in the 
dynamic equation (|2"2")l . As expected, the cusp/groove singularity is completely suppressed by the film-substrate 
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FIG. 1: Surface morphologies of 3% strained films, as obtained from numerical simulations without the wetting effect. The 
simulations start either from small random initial perturbation of a planar film [(a) and (b)] or from a doubly-periodic surface 
profile with wavevector q x o = Qyo = 3/4\/2 and amplitude Ao = 0.01 [(c) and (d)]. Both 3D morphologies, (a) at t = 46.77 
for a portion of a 128 x 128 system and (c) at t = 34.35 for system size A^o x A y o, and also time evolution of 2D cross-section 
profiles are shown. 



wetting interaction, and arrays of strained islands or quantum dots will form and grow. This has been verified by 
our numerical results shown in Fig. [3J where we have used the parameters of e = 3%, h = 0.41, h* nl = 0.3, a w = 2, 
w* = 0.08 or 0.2, and simulation time step At = 1. However, another type of growth instability would occur when 
the strength of wetting interaction is not strong enough (e.g., w* = 0.08 in Fig. [2]), showing as the rapid increase of 
island heights beyond initial time stage and then the blow-up of numerical solution at late times. Such instability 
with unbound growth of island height is absent for stronger wetting effect, such as the effect of w* = 0.2 shown in 
Fig. [5£a) which gives the stabilization and saturation of island evolution. This can be understood from the fact 
that the wetting interaction tends to prevent the depletion of the film-substrate wetting layer and hence suppress the 
mass transport from the valley of an island to its top, leading to the constraint of island height as a result of mass 
conservation. Such effect would increase with the strength of wetting interaction, as can be seen from the results given 
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FIG. 2: Time evolution of 3% strained films with different wetting strength w* = 0.08 and 0.2. (a) Evolution of maximum 
surface height, with a 3D island morphology for w* — 0.08 at t — 330 shown in the inset; (b) 2D cross-section profiles at 
y = L y /2 for to* = 0.08 at t = 330 and w* = 0.2 at t = fOOOO. 



in Fig. [2jb) : thicker film layers between surface islands and shallower valleys are found for larger wetting strength 
w* , as a result of stronger suppression on the valley-to-peak diffusion process. 



VI. RESULTS OF NONLINEAR EVOLUTION 



To examine the detailed evolution of strained film morphology more systematically, we have conducted numerical 
simulations of the full dynamic equation (f2"2"|) for different small film-substrate misfit strains e = 2%, 2.5%, and 3% 
(in such weak strain limit the continuum elasticity approach can be well applied, as shown in most recent studies 4 ^). 
The parameters for the wetting potential are chosen as h*, = 0.3, a w = 2, and w* — 0.2, with all other parameters 
the same as those given in Sec. [V] In our simulations we have used 3 different system sizes L x x L y for each parameter 
set, including the lateral dimensions of 128 x 128, 256 x 256, and 512 x 512, to examine any possible artifacts of finite 
size effects. A numerical grid spacing Ax — Ay — 1 is adopted, and the integration time step is chosen as At = 1. 
The quantitative results given below have been averaged over 20 independent runs for system sizes 128 x 128 and 
256 x 256, and 10 runs for 512 x 512. Also, each simulation starts with a rescaled initial film thickness of ho = 0.67 
for misfit strain e = 2%, ho = 0.5 for e = 2.5%, and ho = 0.41 for e = 3%, all of which are within the corresponding 
instability parameter region for each misfit as determined by Eq. (|25|) . 

Typical simulation results of film evolution and the formation and dynamics of quantum dot arrays are illustrated 
in Fig. [31 for 2.5% mismatch between the film and substrate. At the beginning stage surface undulations occur due to 
film morphological instability determined in Eq. (|25[) . leading to the formation of strained surface islands or quantum 
dots as shown in Fig. [3fa). Note that at different surface locations, islands will form and grow gradually at different 
rates due to the nonlinear effects of elastic interaction. Island coarsening occurs at the next stage, showing as the 
growth of some quantum dots at the expense of other shrinking ones and hence the decrease of island density on the 
film surface. This can be seen more clearly in the corresponding 2D top- view images of Figs. OHd)-(f), which give 
the comparison of island distribution between times t = 1000, 2000, and 10 4 . Such coarsening process becomes much 
slower as time increases, and the system would approach an asymptotic state with steady arrays of strained quantum 
dots. As expected, this late-time state of film surface structures highly depends on the value of film-substrate misfit 
strain, with an increase of island density and a decrease of island spacing for larger misfits. This has been confirmed 
in our results of 2% and 3% films given in Fig. [4] as compared to the 2.5% film shown in Figs. [3jc) and[3]Jf). In our 
simulations no long-range spatial order can be found for quantum dot arrays, even at the late-time stage, agreeing 
with the observation of most experimental and theoretical studies. 

To quantify the above results, we have analyzed the film surface morphology through various time-dependent 
parameters, including the structure factor of surface height, its moments, the maximum height of surface profile, and 
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FIG. 3: Morphological profiles of a 2.5% strained film, at simulation times t = 1000 [(a) and (d)], 2000 [(b) and (e)], and 10000 
[(c) and (f)]. A fraction of 3D morphologies in a 256 x 256 system is given in panels (a)-(c), while the corresponding 2D gray 
scale top- view images of the full system size are shown in (d)-(f). 




FIG. 4: Morphological profiles of (a) 2% and (b) 3% strained films at late time stage of t = 10000. Only a portion of the 
256 x 256 system is shown in the 3D images of (a) and (b). The corresponding 2D top-view images of the full system size are 
given in (c) and (d). 
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FIG. 5: Structure factor of the surface height as a function of wave number q, for different misfits (a) e = 2%, (b) e = 2.5%, 
and (c) e = 3%, system size 256 x 256, and times t = 2000, 6000, and 10000. 




FIG. 6: Time evolution of the three moments of structure factor: (a) mi, (b) mi, and (c) mz, for misfit strains e = 2%, 2.5%, 
and 3%. The solid lines represent the simulation results of grid size 256 x 256, whereas the dashed and dotted lines represent 
the results for grid sizes 128 x 128 and 512 x 512 respectively. Power law fittings (i.e., the thin lines) at the beginning of 
coarsening stage for the 256 x 256 system are also shown. 



the surface roughness, as shown in Figs. [SHU with results of different system sizes presented. The structure factor is 
defined as S(q,t) = (\h(q,t)\ 2 )g, with a circular average over orientation q of the wave number. Typical results of 
S(q,t) are given in Fig. [SJ for different misfit strains of 2%, 2.5%, and 3%. As the misfit increases, the wave number 
related to the peak location of the structure factor becomes larger, corresponding to smaller island spacing and also 
higher quantum dot density as already seen in Figs. H£c),[3jf), and[4jd). Also, for smaller misfit strain larger time is 
needed for the initial formation of islands (if we compare the t = 2000 curves in Figs. [SJa)-(c)), consistent with the 
result of ah in the linear stability analysis. 

Details of the time evolution of quantum dot islands can be characterized by the calculation of various moments of 
the structure factor. The nth moment of S(q,t) is defined as 

Jdqq n S(q,t) 

mnit) = JdqS(q,t) ' (2?) 

which yields the information of the characteristic size and spatial scale of surface structures^. We have calculated the 
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FIG. 7: Time evolution of the maximum value of structure factor, for misfits e — 2%, 2.5%, and 3%. Results for different 
system sizes 256 x 256 and 512 x 512 are shown for comparison. 



first three moments of S(q,t), with time evolution results for three different misfits given in Figs. IH^a) (for mi),[SKb) 
(for 777*2) 1 and[6jc) (for 7773). Three characteristic regimes of film evolution can be identified in each simulation: the 
process of surface instability development and island formation at early times, coarsening of these strained islands at 
intermediate stage, and an asymptotic stage of island saturation at late times. The first two stages can be distinguished 
clearly in the results of moments shown in Fig. [51 as separated by the turning point (i.e., maximum of m„) of the time 
evolution curve for each moment. The increase of m n at the first time stage is due to the continuous appearance of 
new islands at various times as observed in our simulations, and thus the decrease of average island spacing. The time 
range of this early stage of island formation is longer for smaller misfit strain, as a result of overall smaller instability 
growth rate [see Eqs. (|24[) — . This has been verified in Fig. [51 via comparing the three panels of misfits 3%, 2.5%, 
and 2% (from top to bottom) for results of each moment. 

Once most islands have formed (i.e., when the moments m n reach maximum values), they start to coarsen so that 
the average distance between islands increases, leading to the decrease of m n . A power-law behavior of coarsening, 
m n (t) ~ t~^ n , can be obtained, but such behavior is limited to a transient time range at the beginning of the coarsening 
stage. As shown in Fig. [6j this time range is smaller for larger film-substrate misfit strain e, with faster crossover to 
a saturated state. Also, slower coarsening rate has been found for larger misfit, corresponding to smaller coarsening 
exponents /3„ which are identified as (for system size 256 x 256): For e = 2%, ft = 0.1010±0.0008, /3 2 = 0.181±0.002, 
and /3 3 = 0.235 ± 0.003; For e = 2.5%, ft = 0.0702 ± 0.0007, /3 2 = 0.120 ± 0.002, and /3 3 = 0.145 ± 0.004; For e = 3%, 
Pi = 0.0449 ±0.0009, #2 = 0.076 ±0.003, and /3 3 = 0.090 ±0.006. Note that if the structure factor is assumed to obey 
a simple dynamic scaling behavior due to coarsening, one would usually expect that m n (t) ~ t~ nl31 ; i.e., (3 n = n/3i. 
However, the above results of coarsening exponents in the intermediate time range for all different misfit strains do 
not support this assumption, and we cannot identify a simple format of scaling for the structure factor. This might be 
attributed to the complex relaxation of strain energy in the film and the nonlinear elastic interaction between surface 
islands [see e.g., Eq. (1231) ] which are more complicated than that revealed by simple scaling. 

These two regimes of island formation and coarsening are also illustrated in our numerical results for the maximum 
value of the structure factor S max (Fig. [7|, the maximum surface height h max (Fig. [8]), and the surface roughness 
r(t) — ((h — h) 2 ) 1 / 2 (Fig. [5]). The growth of all three quantities can be observed during the first stage of instability 
growth and island formation, which corresponds to the same initial time range as the m n results shown in Fig. [5J 
Both S'max and the roughness r(t) grow exponentially with time at this stage, consistent with the behavior of linear 
instability analyzed in Sec. IIII1 However, at later times during the coarsening process, these quantities show rather 
slow growth and an approach to saturation, even for /i max - Such phenomenon of surface roughness saturation has 
been obtained in a recent study of a nonlinear evolution equation^ 3 -, but not in other studies using different evolution 
equations (which instead observed power-law growth) 34 ' 35 . The limited growth of maximum surface height given in 
Fig. [8] during island coarsening has not been reported in those previous studies, which usually showed a faster growth 
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FIG. 8: Time evolution of maximum surface height for misfits e = 2%, 2.5%, and 3% and system sizes 256 x 256 and 512 x 512. 
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FIG. 9: Time evolution of surface roughness for various misfits e = 2%, 2.5% and 3%. 
different system sizes 256 x 256 and 512 x 512 almost overlap with each other. 



Note that for each misfit, results of 



of /imax such as a power-law behavior—. 

Our simulation results also indicate a crossover from the island coarsening regime to an asymptotic state of steady 
quantum dot arrays, showing as saturated values of S max , h max , and r(t) (see the 2.5% and 3% results in Figs. EH3), 
and more clearly, the saturation of m„ given in Fig. [6] Such crossover can be identified through the slowing of the 
m n decay after the transient of power-law-type coarsening, and occurs earlier for larger misfit strain. To exclude any 
artifacts from finite size effects, we have tested our simulations for different system sizes ranging from 128 x 128 to 
512 x 512, with qualitatively similar results obtained. This phenomenon of the slowing and saturation of coarsening 
process has been found in some experiments of Ge/Si system o 16 ' 17 ' 19 ' 22 and also in some modeling and simulation 
results based on either direct solution of system elasticit y 26 ' 32 or reduced film evolution equation s 33 ' 36 ; however, no 
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FIG. 10: Cross-section profiles of (a) surface chemical potential /i and (b) elastic energy density £, for misfit e = 2.5% and 
different times t = 5000, 7000, 9000, and 10000. The boxed region will be further studied in Fig. [Til 



sign of coarsening termination in annealing films has been shown in some other experimental 20 and theoretical 29,34,35 
work. Also, the effect of different misfit strains on island coarsening and saturation, which is shown important from 
our results here, has not been addressed in most previous studies. Our results given above suggest that much longer 
times are needed to observe the slowing or cessation of island coarsening for smaller misfits, which could be useful for 
addressing the discrepancy of experimental observation: E.g., the suppression or saturation of quantum dot growth 
can be found at relatively short annealing times for Ge/Si(001) system with large misfit (~ 4%) 16 i 17 ' 19 , while it is 
more difficult to observe in SiGc/Si(001) experiments with weak misfit strain (< 1%)22. 



VII. DISCUSSION 



To understand the mechanisms underlying the saturation phenomenon given above, it would be helpful to examine 
the effects of all contributed terms in the evolution equation (f2"2")l . As already discussed in Sec. [V] the wetting 
potential plays a crucial role on the stabilization of surface morphology via limiting the valley-to-peak mass transport 
and thus the growth of island height. Similar mechanism is expected during island coarsening, and the stabilization 
effect of wetting potential should be also important for the island size saturation. On the other hand, the dynamics 
of coarsening involves the redistribution of mass between different islands, a process that cannot be constrained by 
the wetting effect as long as the wetting layer is not depleted. Thus additional factor(s) must be in play to account 
for the island saturation and stabilization process. 

To further illustrate the saturating process, we show in Fig. 1101 the time- varying profiles of the chemical potential 
fx = 7« + W + £ ' at the film surface and also its corresponding elastic contribution (i.e., the surface elastic energy 
density). For simplicity, only the cross-section results at x — L x /2 are presented, for an example of 2.5% misfit film. 
Smaller spatial variations of chemical potential fj, along the film surface are obtained at later times [note the very 
small vertical scale in Fig. [TDT a')]. indicating an approach to an asymptotic saturated state. This is consistent with the 
results given in Figs. [5HS1 for various morphological properties, and also with the evolution profiles of elastic energy 
density £? given in Fig. 110( b). Furthermore, in Fig. 1111 we examine the detailed mechanisms of such saturation 
through identifying the time evolution of various components of chemical potential, including the surface-energy 
contribution 7K and the first and second order elastic energy densities S^f and £ ( 2 ^. We focus on a small region of 
4 islands [see Fig. Illf a)]. representing 3 scenarios of quantum dot evolution: (1) large islands that are growing and 
saturating, (2) small islands that are shrinking [see the middle island in Fig. [TTT a)]. and (3) islands that are migrating 
(see the one at the right corner). As expected from previous analysis^, for an undulated surface (i.e., in the region 
of surface islands), the strain energy is concentrated at surface valleys but released at peaks; the resulting surface 
elastic energy density gradient would drive the diffusion process from the valleys to peaks and thus the growth of 
surface islands. On the other hand, this morphological destabilization process is competed by the stabilization effect 
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FIG. 11: Cross-section profiles of (a) surface morphology and (b)-(d) various energy densities, for 2.5% misfit and times 
t = 1000 (at early stage of instability and island formation), 2000 (island coarsening stage), and 10000 (saturating stage). 
Different components of the film surface chemical potential are shown, including the surface-energy contribution (dot-dashed 
orange curves), lst-order elastic energy density f' 1 '^ (dashed black), 2nd-order elastic density £ ' 2 '^ (green stars), and the total 
elastic contribution £ ! — £^^ +£ (solid brown). 



of surface energy, showing as energy penalty for high-curvature surface areas and hence a spatial distribution opposite 
to that of elastic density (see the dot-dashed curves). This classical view of quantum dot formation has been well 
reproduced in our results of all three evolution stages: the early morphological instability shown in Fig. Illf b) (at 
t = 1000), a coarsening regime in Fig. [TlTc) (at t — 2000), and a saturating stage in Fig. fTlTd) (at t — 10000). 

For islands to be saturated and stabilized, one would expect mechanisms of film evolution involving additional 
energy penalty for large, increasing island size, so that the overall stabilization factors would compensate and suppress 
the destabilization effect (i.e., continuing growth and coarsening of surface islands) caused by stress relaxation. In 
previous studies such factors are usually provided by additional surface energy terms particularly the surface energy 
anisotropy, which has been shown to enhance the surface-energy stabilization effect, constrain the island height, and 
lead to island shape/facet selection and transitio n 16 i 17 ; this effect of surface anisotropy has been deemed essential 
for the existence of steady island arrays in some previous theoretical work (with various assumptions of the wetting 
effect) 26 ' 32 ! 36 . However, in this work we only consider isotropic surface energy. What we identify here is a new factor 
that is due to the contribution of higher-order perturbed elastic energy on the interaction and evolution of surface 
islands, as detailed in Fig. QTJ Positive contribution from the 2nd-order elastic energy density E ( 2 ^ is found for large 
surface islands, showing as an effective energy-penalty term and hence a reduction of strain relaxation effect. [Note 
that this result is still compatible with the well-known strain relaxation mechanism, since the total elastic density E* 
still shows a destabilization effect due to the dominance of lst-order density see Figs. [TTJ (c) and (d).] Such 

effect of E^f becomes important only at late stage with large enough islands, and is negligible for small ones, as seen 
from the comparison between Figs. [TT] (b)-(d). 

To understand this seemingly counterintuitive result which is beyond the conventional view based on linear insta- 
bility analysis, we examine the detailed expression of E^ f which, from Eqs. ([2T])- (P]) . is rewritten as 



E W{r) = 



e* 2 



f(h) + J dv 1 J dv"h{v')G{v - t>, r' - r")h(r") 



(28) 
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where 



f(h) = \Vh\ 2 + u (£ W//e* 2 ) 2 + (1 - V ) (h), 



(29) 




x — Ux(qx - q' x ){q'x + W y ) + q y {q y - q y )(q y + Wx) + 0- - v )q'xq' y [qx(% - q' y ) + q y (q x - q' x 



)]}. (30) 



In Eq. (|28p . the first part /(/i) is always positive, analogous to the "self" elastic energy of a given surface profile that 
serves as a energy penalty to suppress its coarsening; the 2nd part represents the correlation between surface heights 
and thus the elastic interaction between surface islands. Within each island region (particularly near the peak), the 
magnitudes of both parts increase with the island size as verified in our numerical calculations. 

If these 2nd-order elastic contributions are absent or not strong enough, the elastic energy relaxation would increas- 
ingly dominate over the surface-energy stabilization effect, driving the continuing island growth even in the presence of 
the wetting potential. This can be illustrated clearly from our numerical results given in Fig. I12[ where the same film 
evolution equation (|22p is simulated, but with only first-order elastic energy £W/ incorporated. All other parameters 
remain unchanged, including the same wetting potential approximation Eq. The maximum surface height is 

found to increase monotonically with time [see Fig. 1121 (a)], without any slowing or saturation process observed, a 
result that is consistent with previous work^i. Time evolution of the corresponding 2D cross-section surface profiles 
is given in Fig. [121(b), from which two main features of surface dynamics can be identified: (1) Large mass transport 
from film layers to islands is observed, leading to much thinner film layers between surface islands as compared to the 
result shown in Fig. [TT] (a) which incorporates the 2nd-order elastic energy effects. Although the wetting potential 
still has the effect of preserving the wetting layer in-between surface islands and then limiting the diffusion process 
from the depleted wetting layer to the peaks, here such effect becomes relatively weaker as time evolves due to the 
increasing dominance of the destabilization effect of lst-order elastic energy and the absence of "self" energy penalty 
term f(h) for large islands. (2) Mass transport between islands continues to occur, which corresponds to island 
migration or coarsening process and is actually a secondary effect compared to (1). This process cannot be prevented 
by the wetting effect, and can be controlled only by the higher-order elastic energy terms describing island interaction 
and correlation [see Eq. (|2"5)) ]. Thus at late times the island heights increase rapidly, resulting in the formation of 
surface islands with large aspect ratio between height and width as shown in Fig. [T2] The perturbation method used 
here is no longer valid for such high islands, and the simulations will ultimately blow up. This is qualitatively different 
from the results given above with the incorporation of 2nd-order elastic energy, where islands with well constrained 
aspect ratio are obtained which also shows the applicability of the perturbation method developed here. 

All these results indicate that the nonlinearities given by the higher-order strain energy of individual islands and 
the elastic interaction between islands can affect the pathway of film strain relaxation at late evolution times, slow 
down the decrease of total elastic energy via their increasing positive energy contribution for large islands, and 
thus effectively reduce the effect of stress relaxation as the surface instability driving force. Such reduction leads 
to relatively stronger role played by the surface energy and the wetting potential at later times [see the comparison 
between Figs. [TTJ (c) and (d), and between Figs. [TT] (a) and[T3](b)], limiting the mass transport between film layers 
and islands and hence suppressing the island growth and coarsening. 



We have investigated the nonlinear dynamic processes governing the formation, coarsening, and stabilization of 
strained quantum dot islands on the surface of heteroepitaxial films, through the development of a nonlinear evolution 
equation for film morphology. Our study is based on a continuum elasticity model that incorporates the film-substrate 
wetting effect and importantly, on the construction of a perturbation method in Fourier space for determining the 
system elastic properties. In addition to a linear stability analysis which yields the conditions of film morphological 
instability, we have performed large scale numerical calculations of the dynamic equation derived to study the detailed 
behavior of film evolution. We focus on effects of small misfit strains which correspond to relatively large length scale 
of surface nanostructures, and analyze the evolution of strained surface islands/dots using a variety of characteristics 
of film morphology, including the structure factor of surface height, its first three moments, the maximum height of 
surface profile, and the surface roughness. 



VIII. CONCLUSIONS 



16 




FIG. 12: Results of strained film evolution with only first-order elastic energy S^ 1 ^ incorporated. All other parameters are 
the same as those in Figs. 1101 and [TT1 (a) Maximum surface height as a function of time t, with results from calculations up 
to 2nd-order elastic energy (with = £ < - 1 ^ + £^ f ) also shown for comparison; a sample 3D image of island morphology at 
t = 1600 is presented in the inset, (b) 2D cross-section profiles of film surface morphology at different times. 



Consistent with previous experimental and theoretical work, our results have shown three characteristic stages 
of island evolution for post-deposited annealing films, including (1) the early stage of morphological instability and 
island formation, as characterized by the exponential growth of maximum structure factor S^ax and the surface 
roughness as well as the increase of maximum surface height and moments m„; (2) a nonlinear island coarsening 
stage, with a transient power-law behavior of m n decay that appears at the beginning of this stage; and (3) a 
crossover to an asymptotic state of saturated island arrays (although without long-range spatial order), after the 
slowing and suppression process of coarsening. Also, the dependence of these detailed properties on the film-substrate 
misfit strain has been obtained, such as the values of coarsening exponents and the time ranges for the crossover 
between different evolution stages. These have been shown important for the understanding of different, or seemingly 
inconsistent, experimental results particularly for the late time stage of island coarsening or stabilization. On the 
other hand, such dependence does not qualitatively affect our results of the three evolution regimes; same conclusion 
can be drawn for the effect of different finite system sizes used in our simulations. To understand the mechanisms 
underlying the nonlinear evolution of strained films, we have examined the effects of film-substrate wetting potential, 
in particular its role on the suppression of the valley-to-peak mass diffusion process that would lead to wetting layer 
depletion, and its constraining effect on island growth. Furthermore, through a detailed study of time evolution of 
elastic energy density distribution at the film surface, we find that higher-order terms of film elastic energy, which 
incorporate the interaction between strained surface islands and the higher-order "self" elastic energy of individual 
islands, can effectively alter the relaxation pathway of film strain energy at late stage. They play an important role on 
the saturation and stabilization of quantum dot arrays, in particular the crossover to the saturated state with balanced 
multi- island interactions and limited island-layer and between-island mass transport. Thus our results indicate that 
both effects of film-substrate wetting interaction and high-order elastic energy are pivotal for the achieving of steady 
quantum dot arrays and also for the understanding of self-assembly process of strained film heteroepitaxy. 
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Appendix: Third order perturbation results of film elasticity 



As described in Sec. IIII1 our perturbation approach developed here for solving the system elasticity problem can 
be extended to obtain higher order results through a recursive procedure. We have calculated the elastic properties 
of this heteroepitaxial system up to third order, with results presented in this appendix. The 3rd-order perturbed 
elastic energy density is given by 



where 
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Based on the perturbation solutions of the system elasticity, we have obtained the first and second order results for 
clastic stress tensors at the film surface, which are used to calculate the perturbed elastic energy density given above. 
In second order we have 
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while the results for first order perturbation are given by 

^ )/(q) = ~^Bv) {ql + ^ )A(q) ' (A - u) 

^ )/(q) = W^V) {q y + ^ )/l(q) ' (A - 12) 

^(q)=0, (A.13) 

fr«'(q) = 2Ee^(q), (A.14) 

aW(q) = -A i9i A (q)i (A.15) 

aW(d) = -^iq y h( q ). (A.16) 
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